Zeitschrift fiir physikalische Chemie 



Ready, set and no action: A static perspective on potential energy surfaces 
commonly used in gas-surface dynamics 



MSc Vanessa Jane Bukas: Tcchnische Universitat Munchen, Lichtenbcrgstrafie, 4, D-85747 
Garching, Germany 

Tel: +49-89-289-13813, Fax: +49-89-289-13622, E-Mail: vanessa.bukas@ch.tum.dc 



Dr. Jorg Meyer: Tcchnische Universitat Munchen, Lichtenbcrgstrafie, 4, D-85747 Garching, 
Germany 

Tel: +49-89-289-13810, Fax: +49-89-289-13622, E-Mail: joerg.meyer@ch.tum.dc 

Dr. Maite Alducin: Centro dc Fi'sica de Materiales CFM, Ccntro Mixto CSIC-UPV/EHU, 

Paseo Manuel de Lardizabal, 5, 20018 San Sebastian, Spain 

Tel: +34-943-01-8418, Fax: +34-943-01-5800, E-Mail: wapalocm@sq.ehu.es 

Prof. Dr. Karsten Reuter: Tcchnische Universitat Munchen, Lichtenbcrgstrafie, 4, D-85747 
Garching, Germany 

Tel: +49-89-289-13812, Fax: +49-89-289-13622, E-Mail: karstcn.reuter@ch.tum.dc 



Keywords: PES, gas-surface dynamics, dissociative sticking probability, corrugation- 
reducing procedure, neural networks, global minima search 



MS-ID: 

Heft: / () 



joerg.meyer@ch.tum.de 



March 13, 2013 



Abstract 

In honoring the seminal contribution of Henry Eyring and Michael Polanyi who first introduced the 
concept of potential energy surfaces (PESs) to describe chemical reactions in gas-phase [Z. Phys. Chem. 
12, 279-311, (1931)], this work comes to review and assess state-of-the-art approaches towards first- 
principle based modeling in the field of gas-surface dynamics. Within the Born-Oppcnhcimer and frozen 
surface approximations, the O2-Ag(100) interaction energetics are used as a showcase system to accentuate 
the complex landscape exhibited by the PESs employed to describe the impingement of diatomics on metal 
substrates and draw attention to the far-from-trivial task of continuously representing them within all 
six molecular degrees of freedom. To this end, the same set of ab initio reference data obtained within 
Density Functional Theory (DFT) are continuously represented by two different state-of-the-art high- 
dimensional approaches, namely the Corrugation-Reducing Procedure and Neural Networks. Exploiting 
the numerically undemanding nature of the resulting representations, a detailed static evaluation is 
performed on both PESs based on an extensive global minima search. The latter proved particularly 
illuminating in revealing representation deficiencies which affect the dynamical picture yet go otherwise 
unnoticed within the so-called "divide-and-conquer" approach. 
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1 Introduction 



Reflecting the ingenuity of the idea, some 80 years after the seminal con- 
tribution of Eyring and Polanyi [1] it is difficult to envision chemical reac- 
tions without following their footsteps to consider a practical consequence 
of the Born-Oppenheimer approximation (BOA) [2]: potential energy sur- 
faces (PESs). Knowledge of the potential energy as a function of the 
atomic positions provides the basis for a most thorough molecular-level 
understanding of the reaction, in which prominent topological features of 
the PES like minima and saddle points play a central role. Apart from 
such static information, a continuous representation of the PES may also 
be used for explicit dynamical simulations of individual reaction events. 
If allowed by the numerical efficiency of the latter, appropriate averages 
over a sufficiently large number of trajectories can routinely be obtained 
using computational resources commonly available nowadays, thus giv- 
ing access to the kinetics of the reaction and thereby fulfilling one of the 
visions of the pioneers Eyring and Polanyi [1]. 

Although all the more accurate BOA-based numerical solutions of the 
Schrodinger equation have become available during the last eight decades 
[3], the dimensionality of the problem still remains a major limitation to 
employing the PES concept in practice. As demonstrated by Eyring and 
Polanyi [ ] for simple reactions involving triples and quadruples of atoms 
in gas phase, a comprehensive understanding of the PES topology can even 
visually be obtained along two reaction coordinates. For higher dimen- 
sions this becomes an increasingly challenging dilemma which research on 
gas-surface interactions has been battling with for a long time [4, 5]: In 
contrast to gas-phase reactions, the presence of the solid surface destroys 
the overall rotational and translational symmetry of the overlaying ad- 
sorbate species. Consequently, a description based on intramolecular (in 
particular) distances alone is no longer possible [G], thus inflicting an in- 
trinsically higher dimensionality already for diatomics. In the meantime, 
the interaction of localized molecular bonds with the extended electronic 
manifold in particular of catalytically interesting metal surfaces gives rise 
to complex corrugated PES topologies that cannot readily be described 
by a reduced number of reaction coordinates. It is therefore not surpris- 
ing that extensive first-principles based work focusing on prototypical H2 
adsorbates has shown that at least all six molecular degrees of freedom 
need to be explicitly considered in order to arrive at a quantitative un- 
derstanding of central kinetic parameters such as (dissociative) sticking 
probabilities on rigid surfaces [7, 8, 9]. 

An increasing dimensionality is generally also accompanied by an in- 
crease in the number of trajectories required for reliable statistical sam- 
pling, when averaging over trajectories that start from different initial 
orientations and lateral positions of the impinging molecule. The exponen- 
tial growth of computational power has allowed for evaluating thousands 
of trajectories through molecular dynamics (MD) simulations [10, 11] al- 
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ready some 30 years ago - when basing the representation of the high- 
dimensional PES on classical interatomic potentials (CIPs) with "sim- 
ple" analytical forms (like e.g. the Morse potentials used by Eyring and 
Polanyi [ ]). Nowadays, kinetic quantities requiring statistical averaging 
over even orders of magnitude more trajectories can, in fact, readily be 
obtained for (such) numerically convenient PESs, also e.g. as a function 
of initial kinetic energy and varying incidence angles [5] . It soon became 
clear, however, that even the semi-empirical CIPs - with the London- 
Eyring-Polanyi-Sato (LEPS) potential as the most prominent example 
- are severely challenged in capturing the complex bond breaking and 
making involved in a dissociative adsorption or recombination event [11] 
(although promising steps towards more accurate parametrizations from 
first principles have recently been reported [5, 12, 13]). 

On the other hand, ab initio electronic structure calculations, even 
if based on computationally efficient semi-local density-functional theory 
(DFT), have been far too expensive to allow for an on-the-fly evaluation 
of the PES. As heralded by first determinations of sticking coefficients 
through ab initio molecular dynamics (AIMD) [14], this situation might 
just be changing, at least in cases where less extensive statistical sampling 
can be accepted. Nevertheless, the huge number of trajectories required 
to properly determine small reaction probabilities of e.g. < ~ 5% [15] 
practically excludes direct application of AIMD for many systems for some 
time still to come. This is particularly true when taking into account 
nuclear quantum effects, which can be significant for the aforementioned 
prototypical H 2 adsorbates [8, 16]. 

Given this situation, the state-of-the-art is still defined by a divide- 
and-conquer type approach [8, 17, 18]: The 6D PES of a diatomic inter- 
acting with a frozen surface is first mapped discretely by first-principles 
calculations and subsequently "continued in between" (ensuring differen- 
tiability up to at least first order). The actual dynamical simulations 
are finally conducted on this (numerically undemanding) continuous rep- 
resentation by solving the appropriate classical or quantum mechanical 
equations of motion. While conceptually appealing, the practical bottle- 
neck of this approach is the difficulty of procuring and assessing a reliable 
diffcrcntiable interpolation in high (in fact even just six) dimensions that 
usually are very closely intertwined. 

Particularly the latter point shall be the major topic of the present 
contribution. With a number of high-dimensional interpolation tech- 
niques suggested and employed in the context of gas-surface dynamics 
[19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], we will focus here on the 
prevalent corrugation-reducing procedure (CRP) [19, 20, 21] and neural 
networks (NNs) [31, 32] with proper symmetry adaption [26, 27, 28, 30]. 
Providing a nice showcase for the PES complexity often met in gas-surface 
dynamics, we use these two techniques to continuously represent a given 
DFT data set for O2 at Ag(100) from Alducin et al. [ ]. After compar- 
ing and highlighting the methodological differences of both approaches 
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in Sect. 2, we show in Sect. 3 that the two representations yield notably 
different reaction probabilities for the oxygen molecule, despite flawlessly 
passing the typical reliability checks performed to assess the interpolation 
quality. We therefore perform a thorough investigation into the global en- 
ergetics of the two representations both in terms of 2D visualizations and 
of identifying the energetically low lying minima through extensive con- 
figurational sampling. This is found a particularly illuminating approach 
not only in regionally evaluating the accuracy of the PES representa- 
tions, but also in providing an insightful picture into the intricate PES 
topologies. For the O 2 /Ag(100) system it rationalizes the discrepancies 
found for important dynamical observables within the two interpolation 
schemes, thereby highlighting the issue of the interpolation reliability and 
suggesting to generally perform a priori static PES examinations within 
the divide-and-conquer approach. 



2 Methodology 

2.1 Coordinate systems 

As indicated in the introduction and depicted in Fig. 1, the 'full-dimensional' 
representation of the interaction between a diatomic molecule and a rigid 
substrate surface corresponds to a six-dimensional problem. Straightfor- 
wardly, the molecular degrees of freedom can be expressed in terms of the 
Cartesian coordinates of the two adsorbate constituent atoms A and B 

-R cart = { Xa, Ya, A b , Yq,Z b ) , (1) 

Ra. R& 



where for all what is to follow, the origin of the coordinate system is lo- 
cated in the top layer of the Ag(100) surface at the position of a silver 
atom (top site, cf. Figs. 1 and 2). Another representation that more di- 
rectly conveys equivalent configurations due to surface-induced symmetry 
as summarized by Fig. 2 is based on spherical coordinates 



R S ^ = (X,Y,Z, d,0,<p) , (2) 

R 

with the molecular center of mass R, d the molecule's internuclear dis- 
tance, and (■$, if) the polar and azimuth angle, respectively, as defined in 
Fig. I. 

While ij cart can be obtained from R sph without any ambiguities 

Ra,b = R T —j^— d (cos -d sin ip, cos d cos ip, sin -d) , (3) 
we emphasize that special care must be taken with properly defining the 
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Figure 1: Illustration of the 6D coordinate systems used to represent 
diatomic molecules above a frozen surface. Position and orientation 
of the molecule can be equivalently described by Cartesian ij cart = 
(Xa,Y a ,Z a ,X b ,Y b ,Z b ) or spherical R sph = (X,Y, Z,d,i),ip) coordi- 
nates. The right panel depicts the spherical representation of the 
adsorbate-surface system with the axes origin located in the surface plane 
on a top site. The left panel provides a more detailed illustration of the 
molecular orientation tp) within the spherical representation and the 
equivalent Cartesian representation. 
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Figure 2: Examples for symmetry equivalent lateral coordinates on the 
Ag(100) surface given in units of the surface lattice constant a. For dif- 
ferent points (xo,yo) (thick green circle) in the indicated triangular ir- 
reducible wedge (dark gray area), in each panel the equivalents are the 
intersection points (thin green circles) of contour lines of the functions g\ 
(blue) and gi (red) as defined by Eqs. G. Contour values of g± and <?2 are 
given by gi(xo, j/o) an d 32(^0, Ho), respectively. In addition to the transla- 
tion symmetry note the different local point group symmetry around the 
high symmetry sites of the surface: top and hollow sites feature four (left 
panels) , whereas the bridge sites feature only two mirror planes perpendic- 
ular to the surface (right panels) - thus resulting in the aforementioned 
triangular irreducible wedge. Surface atoms are indicated by the large 
light gray circles, with the one in the center of the plot defining the top 
site at which the origin is located. 
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inverse transformation R sph (R c&rt ): 



R 



TOA 
M 



Ra + 



m B 
M 



Rb 



(4a) 



d 



y ftxg ~ X A ) 2 + (Y B - Y A f + (Z B - Z A f 




(4b) 



Z B -Z A \ 180° 



e[0°,180°] (4c) 



arccos 



H(y A - Kb) ' 360° 



+ sign(F B — Y A ) arccos 




180° 



7T 



€ [0°,360°) 



(4d) 



where toa.b is the mass of atoms A, B and M = m A + m B the mass 
of the entire molecule. The distances d\\ and d± are the components of 
the internuclear distance d in the xy plane and along the z axis, respec- 
tively. Particular care must be taken in defining the validity ranges and 
reference directions of the angular variables. The convention that has 
been chosen here (and depicted in Fig. 1) has the configurations with 
•& € (180°, 360°) described by equivalent counterparts in the [0°, 180°] in- 
terval by rotating an additional 180° in ip. Accordingly, the Heaviside (or 
step) function H and the sign function are used in Eq. (4d) in order to 
ensure that configurations in all four quadrants are mapped to the cor- 
rect value of ip € [0°,360°]. Note that permutation symmetry of homo- 
nuclear diatomics is deliberately not included here, thus distinguishing 
R sph (R A ,R B ) ^ R sph (R B ,R A ). The proper one-to-one mapping es- 
tablished here by Eqs. (3) and (4) is crucial for the global minima search 
described in Sect. 2.5 below, as this involves frequent back-and-forth trans- 
formations in order to identify symmetry-equivalent configurations. 

In this context and even more so for the neural network interpola- 
tion described in Sect. 2.3, symmetry adapted coordinates defined (for 
homo-nuclear diatomics) in a similar way to Ref. [30] are of particu- 
lar importance. These map configurations that are equivalent by the 
permutation- and surface-induced symmetries (cf. Fig. 2) to the same 
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Q(i? cart , i? sph ) e 



where 



Oi,3 = 2 S expC-l^.ffLa^.yO 
/e{A,s} 

Q 2 ,4 = n C M-\Z!) ■g ia {X ll Y t ) 
ie{A,B} 

Q 5fi = exp(-^Z) ■ g h2 (X,Y) 
Q 7 = exp(-iZ) 
Qs = d 

Q 9 = [cos(tf)] 2 , 



2tt 

a 

2tt 

a 



2n 



a 

a 



-y 



(5a) 
(5b) 

(5c) 
(5d) 
(5e) 
(5f) 

(6a) 
(6b) 



incorporate the periodicity of the square-shaped surface lattice with the 
lattice constant a. As illustrated in Fig. 2, these have been constructed 
such that symmetry-equivalent points (x,y) in the surface plane result 
in the same pair of function values (g\(x,y), g2{x,y)), with a detailed 
account provided in Ref. [34]. 



2.2 DFT data 

Based on the above spherical coordinate system, a discrete set of ab initio 
total-energy data were calculated within DFT by Alducin et al. in order to 
describe the adiabatic interaction of O2 with a (rigid) Ag(100) surface [33]. 
Briefly summarized, the set comprises data for 12 different so-called elbow 
cuts through the PES, i.e. for molecular configurations of defined lateral 
position of the O2 center of mass (X,Y) and defined molecular orientation 
($,<p). Within every elbow, total energies were calculated over a (Z,d) 
grid, in which the molecule's distance from the surface Z varies from 
to 4.5 A (typically in steps of 0.25 A) for 10 values of the internuclear 
distance d (chosen between 0.9 and 2.5 A). Thus a regular (Z, d) grid was 
constructed containing 2250 energy values of the target PES 

^6D(-R^ ph ) = ^O 2 @Ag(100)(-Ri Ph ) ~ ^Ag(lOO) ~ Eo 2 1 (7) 

and further augmented by values from the molecular O2 binding curve for 
large distances from the surface Z > 4.5 A. The ab initio total energies 
of the clean surface (^A g (ioo))i isolated oxygen molecule (Eo 2 ) and the 
interacting system (i?o 2 @Ag(ioo)) were obtained within a plane- wave ba- 
sis set with a cut-off energy of 515 eV and ultrasoft pseudopotentials as 
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implemented in the VASP code [35, 36, 37, 38, 39, 40]. The exchange- 
correlation energy was calculated within the generalized gradient approx- 
imation (GGA) due to Perdew and Wang (PW91) [41], and the Ag(100) 
surface was modeled in a supercell geometry containing a five-layer slab, 
a (2 x 2) surface unit-cell, and a vacuum distance of 20.89 A [33]. 

2.3 Continuous PES representations 

Based on the DFT data set described above, two different techniques are 
employed within the present work to obtain PES representations that can 
be continuously evaluated along with concomitant forces. The prevalent 
corrugation-reducing procedure (CRP) [ 9, 20] recognizes that a large part 
of the difficulties behind the interpolation of a six-dimensional function 
V(jd describing the PES of a diatomic at a solid surface arises from the 
large energy variations connected e.g. with the strongly repulsive regions 
at short molecule-surface distance. Much of this corrugation, however, 
is equally encountered when the constituent isolated atoms approach the 
surface. The central idea of CRP lies therefore in decomposing the 6D 
molecular PES V^d into a superposition of the 3D atom-surface interac- 
tions 7?-a,b for both atoms (at their coordinates in the diatomic) and a 
smoother function Z CRP , 

V^ p (R sph , R caxt ) = l CRP (R sph ) + K A (R A ) + K b {Rb) ■ (8) 

I p becomes then a more convenient target for the interpolation. In 
practice, this obviously requires additional ab initio data for the con- 
struction of continuous representations of 1Z a,b, which simplifies to a sin- 
gle function in the present case of a homonuclear diatomic. However, the 
lower dimensionality of 1Za,b makes this a much simpler task - including 
a proper incorporation of symmetry [20]. In this like in many preced- 
ing works, Z CRP is then obtained for an arbitrary configuration R sph 
within the range covered by the discrete set {il^ ph } for which DFT ener- 
gies {VQv(Ri Ph )} are available. This is typically done by decoupling the 
six-dimensional problem into four at most two-dimensional independent 
interpolation steps, which can be schematically summarized as follows: 

1. I CRP ({fC h }) (Z ' A):SPli " CS > ^{Z.d- {XuYuKm}) 

2. 2? RP (Z,d; {X^Yu^tpi}) y ' :F °" ricr ) Z 2 CRP (Z,d,^; {X u Y h di}) 

3. I 2 CRP (Z,^; {Xi,YM) ( ^ ):Fb " riCr > I 3 CRP (A,y,Z,d,^; {*<}) 

4. X™ p (X,Y,Z,d,<p; {*,}) l™ p (X,YZ,d,#,<p) 
with X 4 CRP (A, Y, Z, d, 0, <p) = Z CRP (iT ph ). 
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The interpolation functions employed in steps two to four are chosen ac- 
cording to generally expected energy variation in these degrees of freedom, 
and are carefully adapted to its symmetry [19, 20]. 

Neural networks (NNs) as the second prevalent interpolation approach 
are a highly flexible, non-linear model originally inspired by neuroscience 
[32] . Even with utmost mathematical rigor they are shown to be capable 
of approximating any PES to - at least in principle - arbitrary accuracy 
[42, 43]. This motivates to target Vqu directly and equally in its full six 
dimensioniality, allowing to capture intertwinement of degrees of freedom 
which are treated independently in the prevalent interpolation strategy 
hitherto employed in the CRP context described above. Further unlike 
the latter, the input data are not required to fall on regular grids, thus en- 
abling the addition of individual points depending on the desired accuracy 
in certain PES "regions" [28] . Extrapolation capabilities of NNs beyond 
the coordinate ranges of the input data can also supersede those of the 
above described individual interpolations due to the inherent limitations 
of the underlying fixed analytic forms employed in CRP [32]. However, 
this inherent flexibility comes as a mixed blessing: Correct symmetry 
properties can only be guaranteed by presenting symmetry adapted co- 
ordinates to a NN [28, 30, 34], like those introduced in Sect. 2.1. The 
NN thus acts as function J rNN (Q) for the continuous PES representation 
according to 

V^(R cai '\R sph ) = T NN {Q(R cait ,R sph )) . (9) 

In practice and as further detailed in preceding work [27], various multi- 
layer feed-forward NNs of different topologies are fitted (aka "trained" in 
NN lingo) to the DFT input data, after some (166 in the present case) 
randomly selected entries have been moved from the training into a so- 
called test set. During the training iterations (aka "epochs"), the latter 
set is employed to monitor the extrapolation abilities. Using the adaptive 
extended Kalman filter (EKF) algorithm including the modifications de- 
scribed in [27] for training, the sensitivity of the NN is further focused on 
the energetically presumably most relevant parts of the PES by assigning 
training weights to the input data according to 

u Qi (vMR? h )) = aeM-PV 6D (R? h )) ■ (10) 

The parameters a and /3 are chosen in the present work so that 

cj Qi ([-0.25eV;+2.5eV]) = [500; 1] . (11) 

Out of over 70 different attempts we obtained the best fit presented in 
Sect. 3.1 below for a {9-25-25-1 ttZ} NN topology (see e.g. Refs. [27, 28] 
for details on this notation) and the EKF parameters A(0) = 0.9610 and 
A = 0.99670. 
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2.4 Dynamical simulations 

The dynamics on both the CRP and NN PESs are analyzed with quasi- 
classical trajectory calculations that include the initial zero point energy 
of the O2 molecule. A classical microcanonical distribution of the internu- 
clear distance d and its conjugate momentum pd is used for the molecule 
in its ground rovibrational state, which was found equal to 97.5 meV 
[33]. All results are derived from the evaluation of 10,000 trajectories, 
starting at a distance of Z = 9 A from the surface, where the PES value 
for O2 at its equilibrium bond length d eq = 1.24 A is zero in both PES 
representations. 

Along the trajectory calculations, the following reaction events are 
distinguished: 

1. dissociation, when the molecule's internuclear distance reaches the 
value d = 2.4 A (as 2d oq ) with a positive radial velocity; 

2. reflection, when the molecular center reaches the initial starting dis- 
tance of 9 A above the surface and with a positive velocity; and 

3. molecular trapping, when the molecule is neither dissociated or re- 
flected after 15 ps. 

This classification is the basis to arrive at the central kinetic parame- 
ters analyzed in this work, namely the dissociative sticking coefficient and 
molecular trapping probability, defined as averages over the obtained tra- 
jectory data, i.e. the fraction of correspondingly classified trajectories 
over the total simulated trajectories for the given initial kinetic energy. 

2.5 Global minima search 

The low computational cost associated with the evaluation of energies and 
forces on the continuous PES representations allows for extensive config- 
urational sampling in the search for energetically low lying minima. We 
thus employ a rather brute-force scheme that involves a large number of 
individual independent geometry optimizations starting from different ini- 
tial configurations. In these initial configurations the O2 molecule has ran- 
dom lateral positions within the irreducible wedge of the fcc(100) surface 
unit-cell (cf. Fig. 2), random heights and bond lengths between 1.0 — 2.5 A 
and any angular orientation. Each configuration is subject to local geom- 
etry optimization within the Atomic Simulation Environment (ASE) [44] 
until residual forces on the O atoms fall below 0.001 eV/A. The thus opti- 
mized geometries are identified as on-surface molecular adsorption states, 
if the O2 height falls within the range [0, 7.5 A] and its bond length in the 
range [1 A, 2.5 A]. Despite the use of symmetry-reduced initial configu- 
rations, minimization paths can easily lead out of the irreducible wedge 
and thus yield potentially symmetry-equivalent minima. The symmetry- 
adapted coordinates described in Sect. 2.1 were correspondingly used to 
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Epochs 

Figure 3: Evolution of the training and test root mean square error 
(RMSE) as a function of the number of epochs performed during the 
Neural Network training. The inset shows the resulting Neural Network 
energies with respect to the corresponding DFT input data. 

conveniently compare to results from previous iterations and finally es- 
tablish the database of unique adsorption states discussed below. The 
rate with which new such states were added to the set was continuously 
monitored during the sampling run, and used as a criterion to stop the 
search. After 10,000 geometry optimizations on each PES representation, 
this criterion suggested that all inequivalent low-energy minima had been 
identified. 

3 Results &; Discussion 

3.1 Interpolation errors vs. reaction probabilities 

The reliability assessment of a continuous PES representation within the 
divide-and-conquer approach is hitherto commonly restricted to an eval- 
uation of the interpolation or fitting quality. Within the CRP strategy, 
the accuracy is correspondingly checked by comparing interpolated values 
with calculated DFT data not included in the interpolation procedure. 
In the earlier work of Alducin et al. [33] such checks indicated errors to 
fall below 100 meV for a number of different molecular configurations at 
distances Z > 1.5 A. Within the NN approach this kind of performance 
check is already included in the fitting procedure itself by monitoring 
the evolution of the root mean square error (RMSE) with each epoch, 
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cf. Sect. 2.3. For the best fit achieved in the present work the latter is 
depicted in Fig. 3, arriving at train and test set RMSEs of 4.59 meV 
and 8.99 meV, respectively, when terminating the NN training after 150 
epochs. The inset of Fig. 3 shows the resulting deviation of the fitted 
energies to the reference DFT values, which - as expected by the choice 
of the assigned weights - increases in the more repulsive (and allegedly 
dynamically less relevant) regions of the PES. 

The errors obtained with both approaches are en par, if not smaller 
than those reported in a number of preceding divide-and-conquer studies, 
and would generally be considered as reflecting faithful state-of-the-art 
PES representations. As such, the different results obtained in the dy- 
namical simulations on these two PESs and summarized in Fig. 4 come 
rather unexpected: These differences center on the dissociative adsorp- 
tion dynamics of the system, whereas quite similar results are obtained 
for molecular trapping at low incidence energies. For the normal incidence 
case shown in Fig. 4, this specifically amounts to an onset of the dissocia- 
tive sticking probability at about 1.05 eV in the CRP case, while the onset 
is delayed to 1.6 eV for the NN representation. With a monotonously ris- 
ing dissociative sticking probability, this shift in the onset leads to much 
higher values in the CRP case at higher incidence energies, which e.g. with 
0.06 at Ei — 2.0 eV surpasses the equivalent NN probability by more than 
500%. The same trend is observed for various incidence angles tested, 
indicating a much higher reactivity of the CRP representation in gen- 
eral with the dissociative sticking probability already starting to rise at 
significantly lower incidence kinetic energies. 

3.2 Differences in reaction mechanism 

In terms of absolute values for the sticking probability, the differences be- 
tween the two approaches might not appear too worrisome at first sight. 
Due to the large number of trajectories employed in the averaging, they 
are, however, statistically significant. As such, the huge relative discrep- 
ancies question the divide-and-conquer approach in its very core regime 
inaccessible to the alternative direct AIMD ansatz, namely the quantita- 
tive determination of very low reaction probabilities. 

Since both PES representations are based on the exact same DFT 
data, the differences must derive from interpolation deficiencies that go 
unnoticed in the traditional quality indicators based on fitting errors. As 
a first step towards understanding the source of the problem and estab- 
lishing protocols to overcome it, we focus on an analysis of the dissociation 
mechanism. Ideally, this will provide insight into which parts or topolog- 
ical features of the PES are central to the sticking coefficient. In this 
respect, the original work on the CRP PES by Alducin et al. empha- 
sized the role of high energy barriers of about 1.05 eV between bridge 
and hollow sites at a height of 1-2 A above the surface [33]. They re- 
flect molecules with low incidence energies, while molecules with a kinetic 
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Figure 4: Molecular trapping (open symbols) and dissociative sticking 
(filled symbols) probabilities of O2 impinging at normal incidence 0j = 
0° on the Ag(100) surface as a function of the incidence kinetic energy 
Ei. The insets show the position of the molecular center over the unit 
cell, when the dissociating molecules are at Z — 3.5 A (left panels) and 
Z = 1.5 A (right panels) for an incidence energy of 2.0 eV. Note that the 
same coordinate system as in Fig. 2 is used, i.e. surface atoms are located 
at the corners of the depicted unit cell. In all cases, circles and triangles 
correspond to data obtained for the CRP- and NN-PES representation, 
respectively. 
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energy large enough to surmount them are subsequently capable of suffi- 
ciently approaching the surface to finally dissociate in the vicinity of the 
Ag(100) bridge sites. 

Meticulous searches for minimum energy paths performed via the 
Nudged Elastic Band (NEB) method within the present work, revealed 
such barriers to also exist within the NN representation. This suggests 
that the latter yields at least a qualitatively similar reaction mechanism 
as the CRP representation; a suspicion that we find fully confirmed by 
the trajectory data shown in the insets of Fig. 4. Analyzed is the O2 
center of mass position over the surface unit-cell of ultimately dissociat- 
ing molecules as they first reach a height of Z = 3.5 A (left panels) and 
Z = 1.5 A (right panels) for an incidence energy of 2.0 eV. Indeed, in both 
PES representations the molecules predominantly accumulate around the 
bridge position, underscoring the relevance of this PES "region" for the 
dissociation mechanism in both representations. However, in the NN case 
(reflecting the lower dissociative sticking probability) this is a significantly 
smaller number of trajectories. Furthermore, by comparing the distribu- 
tions of the trajectories in both insets in detail, one notices that in the 
CRP case dissociating trajectories are found to be homogeneously dis- 
tributed within the vicinity of the bridge sites, whereas for the NN partic- 
ularly in the direction towards the hollow sites no dissociating trajectories 
can be found. This indicates PES interpolation deficiencies in this area to 
be primarily responsible for the differing reaction probabilities obtained 
with the two representations. 

3.3 Picturing the PES "landscape": a visual compar- 
ison 

A most straightforward way to investigate if and how the two interpo- 
lation approaches give rise to different topologies in those PES regions 
that the preceding analysis identified as most influential for the dynam- 
ical behavior is to visually compare suitable two-dimensional (2D) cuts 
through the two 6D PES representations and include the actual ab initio 
data. Figure 5 depicts three corresponding cuts: A (d-Z) elbow cut over 
the bridge site and for alignment of the molecular axis along the [001] 
direction = ip = 90°), i.e. in the direction of the neighboring hollow 
sites; a ($-<p) angular plot over the bridge site and at a height of 1.5 A, 
i.e. somewhere at the height of the dissociation barriers; and a (X-Y) 
lateral corrugation plot at the same height and for the most favorable an- 
gular orientation identified in the ($-(p) cut (that again points towards the 
hollow sites). Additionally shown are the coordinates of the actual DFT 
data grid in each 2D cut. By construction, the DFT values at these grid 
points are exactly reproduced by the interpolation technique employed in 
the CRP approach (cf. Sect. 2.3). In contrast, within the NN representa- 
tion they are only fitted and the marker size at the grid point reflects the 
corresponding error. 
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Figure 5: Equivalent representations of the CRP (left) and NN (right) po- 
tential energy surfaces. All panels represent two-dimensional cuts through 
the 6D PESs: for X = 0.5 a, Y = and = <j> = 90° (top); for X = 0.5 a, 
Y = 0, Z = 1.5 A and d = 1.75 A (middle); for Z = 1.5, A, d = 1.75 A 
and 6 = <f> = 90° (bottom) . The black points outline the underlying DFT 
grid with the marker shape and size illustrating the corresponding errors 
of the interpolated representation as explained in the figure's legend. 
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Consistent with the weight factors assigned upon training and with 
the initially discussed small global fitting errors, specifically the dynami- 
cally relevant low energy points are very well reproduced by the NN. At 
the high density of DFT data points in the (d-Z) elbow plot, this (or 
the exact reproduction of the DFT values in the CRP case) is enough 
to rather unambiguously determine the continuous representation within 
the corresponding plane, resulting in the remarkable resemblance of both 
plots shown in the upper panels of Fig. 5. This also holds for all other 
eleven elbows contained in the DFT data set (not shown). Despite the 
much sparser supporting grid (of in fact only five inequivalent DFT data 
points), the same still seems to be valid in the angular (#-</?) plane shown 
in the middle panels of Fig. 5. In the lateral (X-Y) plane, however, along 
which DFT sampling is noticeably "poor", both representations exhibit 
significant differences in the levels of corrugation and overall features of 
the PES landscape, with the NN PES featuring an overall more complex 
topology. Similar conclusions are drawn from a large number of 2D cuts 
analogous to those of Fig. 5 for which (some) DFT data points are avail- 
able. These cuts generally reflect a higher degree of "creativity" of the 
NN in describing the PES within areas that are sparsely supported by ab 
initio data. 

In principle, differences in the description of sparsely supported PES 
areas are not surprising, but in fact to be expected: If imagining the con- 
tinuous PES to be like a carpet that has been nailed to the floor at specific 
points (by locally minimizing the RMSE with respect to the DFT data), 
"bumps" and "folds" can freely reveal themselves "between" the (more 
or less) fixed points within both approaches. That there are more such 
"bumps" and "folds" in the NN case is hereby not related to targeting the 
more corrugated Vqd with the NN rather than targeting the smoother in- 
terpolation function I CRP used in CRP. Similar to recent work by Ludwig 
and Vlachos [ ], but including proper treatment of symmetry, we have 
also set an equivalent .F NN ({Q i (.RJ art , i^ ph )}) function as the target for 
the global fitting in symmetry adapted coordinate space with the NN (cf. 
Sect. 2). Pronounced differences in the lateral dimensions, however, are 
still obtained for the different representations. Consequently, the disen- 
tanglement into independent interpolation steps employed in CRP (cf. 
Sect. 2. 3) is instead identified as the cause for the altogether smoother 
CRP-PES. In particular, the Fourier interpolation in (X, Y) does include 
physical assumptions about the PES behavior in these degrees of freedom 
which are absent in the NN case. The aforementioned mixed blessing of 
the latter is thus nicely exemplified: On the one hand, it can produce the 
observed "creative" PES representations. On the other hand, this can 
also allow to better capture more complex PES topologies that are cor- 
relationally induced by all six molecular coordinates and which are then 
beyond the flexibility of the functional form of the Fourier interpolation. 
The latter is e.g. reflected by the lower root mean square error of the NN 
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Figure 6: Potential energy of the local minima configurations detected 
for the CRP (left) and NN (right) PES representations as a function of 
their corresponding distance from the surface, Z. In both cases, differ- 
ent symbols are used to conveniently set apart the low-energy structures 
(crosses) from the more 'shallow' minima located further away from the 
surface (plusses). Note in particular, the NN global minimum (marked in 
black) that has no counterpart in the CRP representation. 



test set described in Sect. 3.1. 

Altogether, the 2D cuts presented in Fig. 5 nicely complement the un- 
derstanding derived from the mechanistic analysis: The (d-Z) elbow cuts 
over the bridge site confirm that both PES representations yield similar 
barriers for this dissociation path that wants the molecular axis parallel to 
the surface and steered into a side-on orientation along the [001] direction. 
In contrast, the markedly different PES topologies encountered in the XY 
plane suggest that the consequences on the sticking probabilities can be 
attributed to the influence of "bumps" and "folds" of both representations 
"in between" these elbows. We turn in the following therefore, to locating 
such topological PES features, and prominently the minima, which as we 
will show not only helps to further pinpoint and understand the source 
of the dynamical discrepancies, but also provides a general protocol to 
assess the PES interpolation beyond the level of fitting errors. 

3.4 PES topology: in search of local minima 

When applying the global search scheme described in Sect. 2.5 to both 
PES representations a significantly different and surprisingly high number 
of local minima is obtained. 21 and 7 bonded (with VjjD < 0) molecular 
configurations are detected for the CRP and NN PESs, respectively, not 
all of which are likely to correspond to physically meaningful adsorptive 
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structures. These minima are plotted in Fig. 6 as a function of their 
vertical distance Z from the surface. We expect minima closest to the 
surface to have the biggest influence on the dissociation dynamics (if there 
is any) , as they should be most effective for an activation of the molecular 
bond, i.e. transfer of energy into the d coordinate. Clearly separated from 
the rest, only one such minimum is observed for the CRP representation 
at a distance of approximately 1.58 A from the surface and an adsorption 
energy of ca. -300 meV. This minimum is scrupulously reproduced within 
the NN representation, which, however, exhibits two additional minima 
in a similar energy range. Since one of these corresponds in fact to even 
the unique global minimum of the NN representation A/nn = (-Xmin — 
0.25 a,Y min = 0.5 a,Z min = 2.UA,d aSa = 1.28A,tf min = 90°,<^ min = 
90°) in the irreducible wedge (cf. Fig. 2), we start by focusing our analysis 
on its influence. 

Returning to the examination of 2D cuts through the PES represen- 
tations, the upper part of Fig. 7 visualizes similarly large differences in 
the lateral degrees of freedom in the vicinity of this NN global minimum 
as discussed already in Sect. 3.3. In order to analyze the influence on 
the dynamics, we have also included (X, Y) coordinates of molecules with 
J5( = 2 eV that pass through the close vicinity of this plane as detailed 
in the figure's caption. Other than in the inset of Fig. 4, we now con- 
centrate on all those trajectories that dissociate on the CRP PES and 
whose equivalent counterparts (i.e. starting from the exact same initial 
conditions) are reflected in the case of the NN. This ensemble of trajec- 
tories (CRP dissociating and NN reflected) does not pass through Mnn 
such that the minimum itself can not directly be held responsible for their 
radically different outcome. Instead, the trajectory distribution is more 
or less centered around repulsive "bumps" only existing in the NN repre- 
sentation shown in Fig. 7 at (Y mi „,0.5 a ± X mia , Z min , d mia , i? min , yw). 
Both "bumps" can be identified with a symmetry equivalent representa- 
tive B NN = (X min ,y min ,Z min ,d min ,z? min ,v3 min - 90°) in the irreducible 
wedge (cf. Fig. 2) with V^, n (Bnn) ~ +0.5 eV. This motivates picturing 
the NN PES representation like in Sect. 3.3 at this "suspicious location" 
within the close vicinity of Mnn also in the other degrees of freedom. 
Fig. 8 compares CRP and NN elbow plots for a corresponding molecular 
configuration. While the CRP PES still shows a dissociation path with 
a barrier of about 1 eV equivalent to the one directly encountered at the 
bridge site (cf. upper part of Fig. 5), its NN counterpart is by more than 
1 eV steeper. Very similar trends are also observed for other elbows corre- 
sponding to different angular orientations of the molecule. For molecules 
therefore with a kinetic energy just sufficient to traverse the minimum 
energy path all these configurations would already lead to dissociation 
within the CRP PES yet be repelled by the much higher barriers on the 
NN PES. 

This fully rationalizes the observed differences in reaction probabili- 
ties shown in Fig. 4 in terms of the differing degree of creativity for the 
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Figure 7: Equivalent two-dimensional cuts of the CRP (left) and NN 
(right) PES representations similar to the lower panel of Fig. 5. The (X- 
Y) plots correspond to Z m i n , d min , •& m - m and ip m i n , i.e. they are cuts along 
the lateral plane through the NN global minimum at (X m j n = 0.25 a, 
Y mia = 0.5 a, Z min - 2.14 A, d min = 1.28 A, i? mln = 90°, p min = 90°) in 
the irreducible wedge (cf. Fig. 2). Note that due to the symmetry of the 
angular coordinates (f? m in, Vmin), a horizontal and vertical mirror plane 
perpendicular to the surface through the hollow site (i.e. the center of the 
panels) is preserved from the clean Ag(100) surface symmetry within these 
cuts. Consequently, there is also a symmetry equivalent counterpart of 
the NN global minimum at (X min = a-X min , Y min , Z min , d min T? min , ip min ). 
Both minima are marked by black crosses. The black circles correspond 
to trajectory data and show the respective positions of the molecular 
center in this cut plane (folded back into the unit cell) for molecules 
that pass closely through this cut plane, i.e. (Z, d, ip) G [i>min>rfmin ± 
0.05 A, tf min ± 30°,^ min ± 20°]. Other than in the inset of Fig. 4), only 
equivalent trajectories starting from the same initial conditions with an 
incidence energy of 2.0 eV are included that dissociate on the CRP (left) 
but are reflected on the NN PESs (right). 
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Figure 8: Equivalent two-dimensional cuts of the CRP (left) and NN 
(right) PES representations similar to the upper panel of Fig. 5. The 
(drZ) elbow plots are shown for X min , Y m - ln , i9 m i n and y m i n — 90°, i.e. 
the coordinates of the "bumps" identified in Fig. 7 backfolded in the 
irreducible wedge (see text). 

CRP and NN pictures discussed previously: While the DFT data grid 
is dense enough to yield the same qualitative dissociation mechanism in 
both representations, additional "bumps" in the NN representation yield 
a much smaller lateral extent of this dominant dissociation pathway and 
thereby lead to a significantly reduced dissociation probability. Whether 
these "bumps" are real or spurious, and correspondingly which PES rep- 
resentation is to be trusted more, cannot be answered without iteratively 
refining the DFT data grid and thus reducing the interpolation freedom. 
However, additional "bumps" and "folds" (aka topological features) nec- 
essarily go hand in hand, and also result in a differing amount of extremal 
points on the PES, as indicated in the present case by the differing num- 
ber of low-energy minima identified on the NN and CRP PESs. Preceding 
any detailed dynamical simulations, a static global minimum search as 
performed here is therefore a suitable and computationally undemanding 
general protocol to assess the interpolation quality beyond the level of 
an arbitrarily chosen test set of ab initio data: Particularly PES regions 
around the identified minima should be supported by a sufficiently dense 
DFT data grid - if not in the original data set, then at least through 
iterative refinement. 

4 Conclusions & Outlook 

In summary, we have presented a systematic comparison of continuous po- 
tential energy surface (PES) representations commonly used in gas-surface 
dynamics within the divide-and-conquer approach. Employing the exact 
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same set of ab initio data [33] for the showcase system O2 on Ag(100) such 
six-dimensional PES representations according to two prevalent schemes, 
namely the Corrugation-Reducing Procedure (CRP) and Neural Networks 
(NNs) , were found to yield largely differing dissociative sticking probabili- 
ties, despite flawlessly passing commonly employed representation quality 
checks. In particular for statistical quantities with low probabilities (such 
as the dissociative sticking probability focused on here) direct ab initio 
molecular dynamics is unlikely to replace the divide-and-conquer approach 
in the near future. This being the case, the reported results are particu- 
larly irritating and dictate the development of protocols to more reliably 
assess and refine the representation quality. 

To this end and complementing an evaluation through explicit dynam- 
ical trajectory simulations we emphasize the value of a static evaluation of 
the obtained PES. In the spirit of Eyring and Polanyi, this concerns fore- 
most the imperative visualization along suitable 2D cuts. In the present 
case this allowed to trace the differences in reaction probabilities back 
to systematic differences of the two approaches in areas that are only 
sparsely sampled by DFT input data. With the latter conventionally 
concentrated in (Z, d) elbow cuts, this applies especially to lateral de- 
grees of freedom, where we found the NN PES representation to exhibit 
a much more complex topology with additional "bumps" and "folds" . In 
contrast, the physical assumptions embodied in the Fourier interpolation 
that is used for these lateral degrees of freedom within the decoupled, 
at most two-dimensional interpolation steps employed in prevalent CRP 
versions necessarily yield a much smoother description. 

For the showcase system, the additional "bumps" were shown to lead 
to a significantly reduced lateral extent of the dominant dissociation path- 
way and thereby to the much reduced dissociation probability obtained 
with the NN PES representation. Both here and in general it is impossible 
to judge which of the two strategies yields the more reliable representation: 
The flexibility of the NN approach, directly targeting the six-dimensional 
PES and its full topology, can allow to foresee/reproduce PES features 
which the CRP approach with its decoupling strategy systematically ex- 
cludes. This can extend to (technical) inconsistencies of the underlying ab 
initio data, which spin-polarized density-functional calculations for sys- 
tems - like the present one - that exhibit a spin transition are particularly 
prone to. On the other hand, this high flexibility brings along the risk 
of a high "creativity" and spurious artifacts in sparsely sampled PES re- 
gions. Ultimately, the only reliable answer can be given by an iterative 
refinement of the input data density in such sparsely sampled PES areas 
that then reduces the freedom of the continuous representation of either 
approach. In practice, the high dimensionality of the problem requires 
focus to be placed on those specific PES areas that are relevant for the 
dynamics. 

In understanding the role of additional "bumps" and "folds" created 
by the NN approach, we find a global minima search to provide a most 
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helpful tool to guiding such a refinement endeavor. In a first place one 
would simply want a continuous PES representation to yield a qualita- 
tively correct topology and therefore knowledge of the representation's 
low-energy minima always comprises a good starting point. Much more, 
however, the added value of locating low-energy regions of the PES lies 
in their general likelihood of influencing the system's dynamics. Even if 
this is not the case - as for the dissociative sticking studied here - the 
aforementioned "bumps" and "folds" (aka topological features) necessar- 
ily go hand in hand. A global search for the latter is in this respect a 
useful (and numerically undemanding) complement to a dynamical tra- 
jectory analysis, which would tediously have to be repeated for different 
incidence energies or angles. The knowledge of the minima can then di- 
rectly be employed to steer a visual PES exploration by focusing on 2D 
cuts that go through their positions. As we have demonstrated here, such 
a look in the vicinity of these minima can reveal "bumps" , which in case 
of sparse DFT data sampling around the latter quickly raises suspicions 
about them being artifacts of the continuous representation. In this man- 
ner, regions of the PES are identified for which adding further data points 
and pinpointing the PES "carpet" will most effectively increase the relia- 
bility of the continuous representation. This holds in particular when not 
having a second such representation at hands for comparison, which after 
all corresponds to the usual working situation. 
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